Magnetized Iron Atmospheres for Neutron Stars 
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ABSTRACT 



. Using a Hartree-Fock formalism, we estimate energy levels and photon cross sections 

, for atomic iron in magnetic fields B ~ 10 13 G. Computing ionization equilibrium 

CN \ and normal mode opacities with these data, we construct LTE neutron star model 



atmospheres at 5.5 < Log(T e ff) < 6.5 and compute emergent spectra. We examine the 
dependence of the emergent spectra on T e g and B. We also show the spectral variation 
with the angle between the magnetic field and the atmosphere normal and describe the 



significant limb darkening in the X-ray band. These results are compared with recent 
detailed computations of neutron star H model atmospheres in high fields and with 
low field Fe and H model atmospheres constructed from detailed opacities. The large 
spectral differences for different surface compositions may be discernible with present 
X-ray data; we also note improvements needed to allow comparison of Fe models with 
high quality spectra. 



Introduction 



Neutron star surface emission has long been sought as a measure of the star's thermal 
history and, hence, a probe of the equation of state of matter at supra-nuclear densities. Since 
temperatures remain near 10 6 K for ~ 10 6 y, this emission is best probed in the soft X-ray band. 
Recent detections of a number of neutron stars with ROSAT (Becker 1995) show good promise of 
enabling a study of neutron star cooling. While spectra are still generally fit with blackbodies at 
present, it has been shown (Romani 1987) that radiation transfer through the surface layers can 
have very significant effects on the emergent spectrum, depending on the composition. 

Conditions on the neutron star surface are very uncertain. While the supernova explosion 
has a mass cut in the iron layer, envelope fallback, accretion during neutron star evolution, and 
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possibly spallation by energetic magnetospheric particles can introduce light elements (eg. H and 
He) to the neutron star surface. Because gravitational sedimentation is rapid, the lightest element 
present should dominate the photosphere and the emergent spectrum. This uncertainty makes it 
important to compute emergent spectra for both light element and heavy element atmospheres. 
As shown by Romani (1987) the general trend is for light elements, with Kramer's law opacities 
at high E, to have emergent spectra with large excesses in the Wien tail over the equivalent 
blackbody. Heavy element atmospheres give spectra much closer to blackbodies, but show strong 
line and edge features. Recent updates of low B model atmosphere computations (Rajagopal & 
Romani 1996 (RR96); Zavlin, Pavlov & Shibanov 1997 (ZPS96)) provide a good set of models for 
comparison with low field neutron stars. As an example, RR96 found that H spectra fit ROSAT 
PSPC observations of the low-B millisecond pulsar J0437-4715 much better than did heavy 
element atmosphere spectra. 

Most observed neutron stars, however, have surface fields at least as strong as the surface 
dipole, 10 125 G or higher. Detailed treatment of the atomic cross sections and the equation of 
state has allowed Pavlov and colleagues (Pavlov et al. 1995; Shibanov et al. 1992) to study 
emergent spectra from magnetized H atmosphere neutron stars. Unlike the recycled pulsars, 
though, the young, high field neutron stars have presumably experienced little accretion. It is 
therefore important to estimate the spectral effects of heavy element atmospheres to compare with 
the magnetic H results. A first attempt was pursued by Miller (1992) who obtained approximate 
energy levels and cross sections for heavy atoms in strong fields (Miller & Neuhauser 1991) and 
computed emergent spectra. These opacities were very approximate, employing only bound-free 
cross sections and averaging the polarization modes. Pavlov et al. (1995) argue that these 
approximations are too severe for the high field models. 

In this paper we present substantially improved magnetic Fe atmospheres, incorporating 
more detailed opacity and equation of state estimates, normal mode radiative transfer and more 
realistic magnetic geometries. These approximate atmospheres will provide a useful baseline for 
comparison with the magnetic H results. Since fully detailed Fe atmospheres will be very difficult 
to construct, we focus on the broad spectral differences between models and note the areas where 
more detailed models can improve the spectral estimates. 

2. Atoms in a High Magnetic Field 

When the electron cyclotron radius p = 2.5 x 10 _10 (S/10 12 G) -1 / 2 cm becomes smaller 
than the approximate zero- field outermost electron orbital sizes of an ion a ~ ao(Z c g + l) -1 , 
(ao = 5.2 x 10~ 9 cm, Z e s the net charge of the ion), the ion's atomic structure perpendicular to 
the field is altered significantly. For hydrogenic atoms of nuclear charge Z e g + 1 in the ground 
state, this occurs at a magnetic field B = (Z c g + l) 2 B c , where B c = 2.35 x 10 9 G. Well below this 
field the electron wave functions are approximately spherical, whereas far above it they are best 
represented as cylinders given by ip mv (z, p,(f>) = fnmu(z)& n m(p, 4>), where (p,4>) are cylindrical 
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co-ordinates about the magnetic field direction z, Q nm (p,4>) are the cylindrical Landau states, 
and fnmv{ z ) is the wave function in the direction of the field. This assumption of separability in 
cylindrical co-ordinates is sometimes called the 'adiabatic approximation'. Here u and m are the z 
and <p quantum numbers, and n is the Landau level. The <&o m (p,4>) wave functions are annular, 
with radius a± = y/2m + 1/5 and width ~ p (e.g Meszaros 1992). 

To compute energies and electron configurations for iron in strong magnetic fields, we use 
a multiconfigurational Hartree-Fock code developed by Neuhauser, Langanke, & Koonin (1986; 
see also Miller & Neuhauser 1991). This code computes the wave functions and state energies 
in the 'adiabatic approximation', with Landau level n = only, since the first Landau energy 
(cyclotron energy, E c = 11.6-Bi2keV for a magnetic field of 10 12 -E>i2G) greatly exceeds both kT 
and all the photon energies E~ we consider. For the innermost electron of iron at the magnetic 
fields we consider, B = 10 125 G and B = 10 13 G, this code gives binding energies that can be 
low by as much as 15% because the Coulomb field of the nucleus is too strong for the separability 
assumption. However, for electrons with binding energies less than ~3 keV, whose transitions 
are most important for the observed spectrum, the code gives individual orbital energies accurate 
to better than 1%. For each ionization state, we have computed the energies and electron 
configurations of the ground state, and the first several (~ 10) excited states which generally 
extend to > 5kT for the typical local atmosphere conditions. As electrons are added to the bare 
nucleus, in the ground state they occupy the "tightly-bound" u = states, filling m = 0, 1,2, ... 
in order. These states have no 5 node, and hence are localized nearer the nucleus than the 
"loosely-bound" v > states. Ground states with many electrons also fill the v = 1, m = 0, 1 
orbitals, with only one electron per orbital (a spin flip requires 5E ~ E c ). The excited states we 
have are all formed by promoting the outermost electrons to higher m, which costs little energy 
because of the weak m dependence of a± (above). We use these configurations to calculate the 



bound-free absorption cross section in each of the (+, — , z) polarization modes (described in §4.1) 
as a function of frequency for each individual electron, for at least the lowest five energy states of 
each ion. We also calculate the energies of the most important configurations accessible from the 
ground state via dipole allowed transitions, and some characteristic oscillator strengths for these, 
for transition energies up to 10 keV. These configurations are made by promoting any one of the 
outermost (large-m) electrons to v > 0. 



3. Ionization Equilibrium and the Equation of State 

Let the r^-ionized atom (Z e g = r) have excited states k with energy levels E r ^. The 
ionization energies are Xr = -EV-,0 — ^(r-l),0) while the excitation energies are e r ^ = E r ^ — E r $. 
Then ionization equilibrium in a strong magnetic field gives the ion densities n r as 

„ _o ra (r~i) f m e kT \ 3 / 2 Vr sinhr^.!) Ve ^_ Yr/fcT E^o e~^ fc/fcT m 
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2-kK 2 J sinhr/ r rj(r-i) tanhr/ e J2T=o e e( - r -^' k ^ kT 
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(Khersonskii 1987), where the sums are over the excitation states of the individual ions. Here 
rj e = hujB/2kT and r/ r = hto r /2kT, where wg = eB/m e c is the electron cyclotron frequency, 
co r = reB /mp e c is the cyclotron frequency of the r'^-ionized atom, m e is the mass of the electron, 
and wiFe is the mass of an iron nucleus. The sums are the partition functions for each ion, 
and the relative abundance of each excitation state is given simply by its Boltzmann factor 
divided by that sum. At high densities in our atmospheres, the electron continuum energy can be 
significantly depressed; we approximate this depression by decreasing the ionization energies by 
E dep = (3/5)e 2 (4vrp/3m Fe ) 1 /3((Z e 2 fr )/(Z cfr )) = ll.Sp 1 / 3 ((Z^ a ) / (Z cS )) eV, with p the atmosphere 
density in g/cm 3 (Cox & Giuli 1968). Since we do not have an exhaustive list of e r ^ to compute the 
partition function sums in (|l|), we take a constant level spacing from the highest available excited 
states to complete the partition sum as a geometric series. This method does not include the 
loosely-bound states v = 1,2,3..., of which we estimate there are ~ 10 for each m, spread between 
~ 1/3 the tight binding energy and the suppressed continuum. At the highest temperatures 
present in our atmospheres, the Boltzmann suppression for these states only partly compensates 
for their large number, and they may contribute significantly to the partition function. Under 
these conditions, more detailed treatments of ionization equilibrium will require exhaustive line 
lists, and a means of calculating pressure effects on partition sums and ionization tailored to the 
magnetic plasma. 

Our equation of state under these extreme conditions is of necessity approximate, despite 
some justified simplifying assumptions. Although the gas deep in our atmospheres would be 
degenerate at low fields, with high B the degeneracy is strongly suppressed. In particular, the 
gas is non-degenerate if T 3> Tp = 0.6(p(Z c g /26)) 2 / B\ 2 (Hernquist 1984). For our atmospheres, 
B\ 2 = 10 — 100 and p(Z e g /26) < 10 3 , so the gas is far from degeneracy throughout (Fig. 2). 
Also kT <C ItjjOb] in the adiabatic limit, this ensures that only one Landau level is occupied, 
although more exact treatments require the consideration of higher n (Potehkin, Pavlov and 
Ventura 1997). The plasma parameter T = Z 2 ff e 2 (47rnj/3) 1 / 3 /fcT, where m is the ion number 
density, is far below the solidification value (r ~ 170). While our atmospheres do extend well 
past r = 1 at high T e ff, we ignore the increasing non-ideality under these conditions (eg. Fehr 
and Kraft 1995). Accordingly, we solve for the ionization equilibrium conditions iteratively, using 
(|l]), to find the number density of all species (and hence p(P,T)) required by the ideal gas law 
to produce the pressure P(t) needed for hydrostatic equilibrium (§||). The mean ionization level 
(Z e ff) is shown for selected r in Fig. 2. We have tested the sensitivity of our final results to the 
omission of non-ideality corrections to the equation of state in two ways. First, we have adopted 
20% of the Debye-Hiickel correction as approximately appropriate for the densities near the base 
of our atmospheres (Rogers 1981), to estimate the structure perturbation in our 10 6 K atmosphere 
model. Holding the temperature run fixed, we find that although there are ~ 5% variations in 
density at the largest depths considered, at the surface the differences in the emergent spectrum 
are below ~ 1% near the peak, and the total flux difference is ~ 0.25%. We have also made the 
extreme assumption that the plasma becomes an incompressible liquid at T = 3 and tested the 
effect on the emergent spectrum: ~ 5% changes were seen in the emergent flux near the weakly 
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bound lines, but significant continuum variation was found only for EL > 20kT c g. We therefore 
ignore plasma corrections in the following computations. 



4. Opacities 

4.1. Radiation Propagation in the Magnetized Plasma 

In a magnetized plasma, radiation may propagate only in two specific polarization modes, 
called the normal modes. Their properties are obtained in general by finding the polarizability 
tensor of the plasma, and using it to solve Maxwell's equations for the case of a plane wave. For 
a collisionless magnetized plasma, the resulting mode polarization components are easily given 
(Ginzburg 1970) in the co-ordinate system (x',y',z'), where z' is in the propagation direction k, 
and the magnetic field direction B is in the y'-z' plane, at angle 9 from the z' axis. The relative 
strengths of the components for modes j = 1,2 are 



2u 1 / 2 (1 - v) cos/ 



e j , usm 2 9 + {-iy- 1 {u 2 sm 4 9 + 4u{l-v) 2 cos 2 9) 1 / 2 . , 

1/9 W 

j . u I vsm.9 w cos 6* sin # 

g"' ^ — ^ g-^ j g"-* 

z u — (1 — v ) — uv cos 2 9 x u — (1 — v) — uv cos 2 9 y ' 

where u = (wg/w) 2 , and v = {uj p /uj) 2 with u> p = 4irn e e 2 /m the plasma frequency. 

We normalize these components, and rotate them to the basis (x,y,z) with z in the B 
direction, and k in the y-z plane, retaining the definition of 9, not to be confused with the angle 
9 n between k and the atmosphere's outward normal. Finally we construct e J + = (e? x + ie J y )/\/2 and 
e?_ = [e? x — ie{i)/y/2, the polarization components in the rotating co-ordinate basis in which the 
opacities are most easily calculated. The indices of refraction are 

n 2 = 1 2 < l ~ V -l (3) 

3 2(l-v)-usm 2 9 + (-iy(u 2 sm' i 9 + 4u(l-v) 2 cos 2 9) 1 / 2 K ' 

For the large u we consider, the j = 2 mode has an imaginary index of refraction (n = in, n real) 
below the plasma frequency, when 1 < v < sec 2 9. In this case we add the extra "skin-depth" 
opacity Tg(u) = noj/c to the total mode opacity of §[0|. In practice all but a negligible fraction of 
the emergent flux from our atmospheres is above the plasma frequency at the depth of spectrum 
formation. 

In fact, absorptive terms in the polarizability tensor corresponding to the various opacities 
affect the mode decomposition. As a test, we treat electron-ion collisions as in Ginzburg (1970), 
approximating non-magnetic free-free absorption: the effect is negligible except well below the 
plasma frequency. Nonetheless we expect that proper inclusion of the atomic opacities will change 
the mode polarizations somewhat, especially in the line regions. To date, the only treatment of 
this effect is for hydrogen (Bulik &: Pavlov 1996). 
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4.2. Thomson Scattering and Free-Free Absorption Opacities 

Meszaros (1992) gives the normal mode electron scattering and free-free opacities as 



a = (<TT, Off) 



UJ 2 , ,, UJ 2 



e+\ 2 g± + -7— r^\e-\ 2 g_L + \e z \ 2 g\\ 



[lo + lob) 2 (w-Wfl) 



(4) 



For Thomson scattering, ax = (8ir /3)ap(h/mc) 2 is just the non-magnetic cross-section, where op 
is the fine structure constant, and g = 1. For free- free absorption (Brehmsstrahlung), 

A 2 3^ c2 /7 2\ l-exp(-faj/fcT) 
a// = 4^ a£ (Z eff )n;^- 1/2 , (5) 

is also the field-free cross-section, and the Gaunt factors are given by 

/ rp r n 1 P r 2// ,^ 1 Cl(a+) + gl(Q-) 

y|,( W ,r,B) = J_jM-P /(2mkT)] {p2 + 2mhw)1/2 (6) 

a± = (p ± [p 2 + 2mhw] 1/2 ) 2 {2mhuj c y 1 
Co (a) = exp(a)-E 2 (a)/a Ci(a) = exp(o)[£?i(a) — -E 2 («)] , 

(e.g. Meszaros 1992), where E\ and i?2 are the exponential integral functions. Under the 
conditions of our atmospheres, these Gaunt factors are ~1 to 10; they slightly steepen the inverse 
dependence of opacity on frequency. 



4.3. Bound-Free and Bound-Bound Opacities; Line Broadening 

Taking the bound- free (+, — , z) basis cross sections from the Hartree-Fock wave function 
sums, we add the contributions from all electrons in each ion. The result for any excitation state 
of an ion is a bound-free edge and a falloff with complex structure. In the case of hydrogen, 
these features can be affected by the thermal motion of the ion transverse to the field (Kopidakis, 
Ventura & Herold, 1996); and, if adiabaticity is not assumed, by resonances with states with 
non-zero Landau level (Potekhin, Pavlov & Ventura, 1997). We expect the former effect to be 
much diminished for iron, while the latter could be as significant as other inaccuracies inherent in 
the 'adiabatic approximation' (see §§). When atoms are excited to a state for which an explicit 
bound-free cross section is not available, we assign them the cross section of the highest available 
state. 

To estimate the line opacity we use a list of energy levels from the Hartree-Fock modeling. 
Because of uncertainty in the energy level computations, line energies may be in error by as 
much as 10%, but the general distribution of bound-bound transitions is representative of the 
true spectrum. We focus here on dipole allowed transitions, which obey the selection rules 



-7- 



Am = ±1 for even Av, and Am = for odd Av (Ruder et al. 1994). We find that in the 
adiabatic approximation, the energy level spectra for individual ions are roughly self-similar. 
Therefore, for each ionization state we use the Am = 0, Av = 1 transition of the outermost 
ground state electron, which we always have, to (linearly) re-scale the most extensive line energy 
lists available (for Z e R = 2 at the lower field, Z e Q = at the higher) and apply it to the ion. 
For most ions we have explicitly computed a few allowed transition energies; these are always 
used, but in any case they match well to the predicted energy from the re-scaled line lists. 
We have attempted to compute a 'complete' set of dipole allowed transitions from the ground 
states only; we generally do not have explicit energies for all relevant bb transitions from excited 
states. To approximate the total distribution of oscillator strength from these transitions, we 
estimate the number of missing lines for each excited state and assign the ground-state line 
spectrum to the unaccounted-for fraction. While we have not computed the oscillator strengths 
for all transitions in detail, we can estimate the line opacity by noting that the strengths, 
defined by / a(E)dE = (2ir 2 e 2 h/mc)f osc , take on characteristic sizes for the different transitions. 
Furthermore, the type of transition determines to which component of the polarization basis the 
line couples. We use / osc (couples to + mode) 0.06 (Am = 1, Av = 0); f OS c(z mode) ps 0.8,0.06 
(Am = 0, Av = 1,3); and / osc (+,-) ~ 5 x 10~ 4 ,5 x 10~ 5 (Am = 1,-1; Av = 2) for the 
oscillator strengths at B = 10 12,5 G. For B = 10 13 G, the corresponding strengths are / OS c(+) ~ 0.04 
(Am = 1, Av = 0); f osc (z) ps 0.4, 0.03 (Am = 0, Av = 1, 3); and /«(+, -)w9x 10~ 5 , 9 x 10~ 6 
(Am = 1, —1; Av = 2). These oscillator strengths are accurate to within a factor of ~ 2, which is 
sufficient for our modeling needs. 

We include two very different broadening mechanisms in our calculation. First, at the high 
densities and temperatures characteristic of neutron star atmospheres, line features experience 
significant collisional broadening. For an electron in an orbital of characteristic radius a about 
a nucleus of charge Z, the energy shift due to the dipole induced by an external electric field F 
is AE ~ F 2 a 3 /Z. For collisional impact, we have F = e/r; 2 mpact , and Z = Zy mc the net charge 
interior to the electron absorbing the photon. We use the Weisskopf method to approximate the 
broadening, since the impinging electrons are confined by the magnetic field to move in a straight 
line, giving a Lorentzian cross-section with width 



^r coll = vr 5 / 3 2V3^/3 a 2 ^3(^)1/6^ ^ 41 i/eK/HM 

m e 2j J Z ' 

line line 



i 2 



{a /Zi 



line j 



eV , (7) 



which comes from the amplitude and r-dependence of AE (see Mihalas 1978). For a we use the 
largest dimension of initial or final state, with dimension a± transverse to the field as given in §||, 
and a z parallel to the field estimated as follows. Orbitals with v = are localized near the nucleus; 
we need their a z only if it is bigger than a±, in which case it is given, through minimizing the 
total energy of an annular charge distribution, by a z = ao/ ln(ao/aj_) (e.g. Meszaros 1992 §2.4). 
For v > we estimate a z = (ao/Zn ne )[(v + l)/2] 2 , which comes analytically from the Schrodinger 
equation for fo mv (§§) for odd v states of hydrogen at high field (Ruder et al. 1994, §6.4.4). Since 
energies vary smoothly with v regardless of parity, both in that work and in ours, we use the 
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above estimate for even v also. 

Secondly, the thermal motion of the ions past the large B also causes significant energy level 

perturbations (Pavlov & Meszaros 1993). We can approximate this effect by a Lorentz electric 

field F = (v±/c)B, where v± is the velocity transverse to the magnetic field. As in the case of 

hydrogen, we take this 'magnetic broadening' to be one-sided, since the Lorentz field always results 

in elongation of the otherwise circular transverse wave function. We estimate the energy shift as 

above by AE ~ F 2 a 3 /Z\- me , taking a to be the larger a± of initial and final state, since F must be 

7] 



transverse to B. Using vj_ ~ ^j2kT/m-p e we get F = 5.8 x 10 Bi2\JTq cgs, yielding a width of 

hT mag = 0.035(2m + if' 2 I^phl e v . (8) 

•^line 

Making proper allowance for the higher thermal velocity of hydrogen, equation (|J) agrees with 
the numerical example in Pavlov & Meszaros (1993) within a factor of 2. Note that since 
it is independent of density, magnetic broadening is most important near the surface of our 
atmospheres. We add the two broadening mechanisms, giving a redward half-Lorentzian with 
width KTn = hT ma g + hT co \i, and a blue line wing with width HT^ = KT C o\\- The line profile, 
normalized to the / osc as above, is thus 

^-^{^( wm^? )^ E< - E ° (9) 

where Eq is the line energy. As this net broadening is sufficient to blend many of the lines, the 
exact values of the line energies are not crucial in estimating the b-b contribution to the overall 
opacity. 



4.4. Total Radiative Opacity 

When our atmosphere calculation calls for the radiative opacity in each mode at a given 
angle and frequency, for some ambient temperature and density, we sum the scattering, free-free, 
bound-free and bound-bound (+, — , z) cross sections for that frequency over all ionization 
and excitation states weighted by ([[]) at the appropriate temperature and density, giving 
<T +i _ t z (u, 9, T, p). The mode cross-sections are then given by 

aj(u,9,T,p)= J2 4(^0) 2 °i(",T,p) j = 1,2 (11) 

i=+ , — , z 

We convert these from cm 2 /particle to cm 2 /g, and multiply by the density to get the opacity 
Kj(u>,0,T, p) in cm -1 . 
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Fig. 1. — Left: Contributions to the radiative opacity at 8 = 55°, at total mean optical depth 
(defined in §^) tt = 1 in our B = 10 12 ' 5 G, T e g- = 10 6 K, &b = 0° Fe atmosphere. The two modes 
are shown by solid and dashed lines respectively. Note that significant b-f opacity is present at 
large E. b-b transitions between tightly bound states appear at low energy; tight-weak transitions 
are at higher energy. Right: Same for B = f0 13 G. 

To illustrate the relative importance of different contributions to the opacity, we show 
opacity spectra for the two normal modes at 55° to the magnetic field, under characteristic 
atmosphere conditions. Free-free opacity dominates at low energies, but b-f and b-b processes 
contribute significant opacity at higher energies; electron scattering is never significant in these Fe 
atmospheres. 



The magnetic field, at an angle Qb to the atmosphere normal, causes significant anisotropies 
in the energy transfer. In addition, ZPS96 have shown that even without a field there can be 
substantial limb darkening and that the emergent spectrum will vary with /i = cos 6 n the cosine of 
the viewing angle to the atmosphere normal. For emission from a small region on the star, these 
effects will cause substantial variation with the viewing angle. We therefore follow the angular 
dependence of both internal radiation field and emergent spectrum, although our atmosphere 



5. Energy Transfer in the Atmosphere 



5.1. 



Radiative Transfer 
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temperature structure remains plane parallel. 

We define a grid of angular co-ordinates (6,(j)) about the magnetic field direction. Given 
the state of the atmosphere at each depth, we use the opacities of §3 to compute frequency- 
and 9- dependent optical depths, by summing Tj(d,uj,9) = J Kj(u),9,T, p)dz from the surface to 
atmosphere layer number d, for each theta in our grid. This precedes any consideration of the 
atmosphere normal direction, and the sums are done using the normal thickness dz of each layer. 

We next compute radiation fields at each 9 and 4> at a given atmosphere level. The intensity 
of outward going radiation is given by 

roc 

I j [ Tj (u,e),n] = [BUr')/2}e-^-^^dT'/f, (0 < M < 1) (12) 
and inward going radiation by 

IjfcfaO),!*] = [ Tj [BUr')/2]e-^-^-^dr'/(-^ (-1 < M < 0) (13) 
Jo 

where \x = cos 9 cos 0_b + sin 9 sin b cos (f> (taking the normal direction to have <j) = by definition) 
and B W (T) = (^ 3 /[ 27I " 2 c 2 (e^ /A:T - 1)] is the Planck function, half of which gives the source 
function in each mode for LTE. When Arj(uj,9) > 20 in a given layer we use the diffusion 
approximation for that frequency and mode: I[Tj(uj, 9), fj] = B u [Tj(u,9)] + [idB u /dTj{oj,9). The 
intensities are then summed to find the astrophysical flux Fj w = (l/n) J J fj,Ij(u,9,(f>)dO,. We find 
that 18 zones in 9 and 9 zones in <ft (zero to it only, by symmetry) give better than 0.1% accuracy 
in the intensity integral even for test functions with rapid variations in 9 and <j>. To monitor 
atmosphere flux the above computations are performed on 200 energy bins, logarithmically spaced 
between leV and lOkeV. 

5.2. Electron Thermal Conductivity 

Electron thermal conduction provides a significant component of the heat flux at large depth 
in our coolest atmosphere models. The electron thermal conductivity is significantly modified 
under high field conditions (Hernquist 1985). Because the electron gas in our atmospheres 
is non-degenerate (see §||), and occupies only the lowest Landau level, the conductivity Ay 
(Hernquist's Kn) along z can be readily computed by evaluating the integrals in Hernquist's Eq 
(3.5), using (3.8) and (3.10). The chemical potential is given by (2.12), and its derivative required 
by (3.10c) is taken analytically. Electron conduction perpendicular to B will be much smaller, 
and we take it to be zero. We find that in our coldest B = 10 13 G atmosphere, in which it is 
most significant, the electron conductivity carries up to 80% of the heat flux at T ra d = 10 (defined 
below). The extra mode of transport affects the temperature structures and emergent spectra 
significantly. 
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6. Converging the Atmosphere and Obtaining Spectra 

Deep within the atmosphere, the optical depths are large for all frequency bins, and radiation 
flows in the diffusion limit with flux = (1/vr) JJ /j,(9, (j)) 2 (&B W / dr^e) svaOdddcj). We therefore 
define a radiative mean opacity 



7T 



4aT 3 




1 



(14) 



equivalent to the Rosseland mean opacity; the corresponding optical depth scale is rrad. The 
total mean opacity is then given by k^} = K~^ d + where k c = 16oT 3 /(3A||Cos0b) gives the 
effective conductive contribution. This total opacity defines a useful depth scale tt < Trad, as the 
temperature structure tends to the grey solution at large depth on this scale. Our atmospheres 
have 146 depth zones with —6 < log tt < 2; the very low r layers are required by the extreme 
frequency and angle dependence of our opacities. 

As usual we relate the optical depth zones to physical depths by summing the equation of 
hydrostatic equilibrium down from the surface: dP = [g s p/KT(Tr)]dTT. Because kt = kt(p,T), 
we solve iteratively for the physical zone size (c/. RR96). We then evaluate the total radiative 
flux F ra d = J2j=i.2 I Fjujduj as above at each depth zone, and add the conductive heat flow 
F c = (A|j/7r)dT/dz to obtain the total energy flux at each depth. This is adjusted to the constant 
flux solution, using the Lucy-Unsold temperature correction scheme, as in RR96. We iterate until 
flux errors are < 1% throughout the atmospheres. 

In the surface layers, we also consider the behavior of / = J du> j sin9 d#d</>X^=i K j(I ~ B w /2). 
Although convergence to / = is a useful diagnostic, we find that a simple A iteration scheme 
using /, as suggested by Shibanov et al. (1995), provides no help in reaching convergence. Full 
scale Feautrier schemes could, however, speed convergence considerably. 

As noted above, both angle dependent specific angular fluxes and the total emergent spectrum 
are of interest. These are computed from equation (12) at the surface layer at 10 3 log-spaced 
energies for higher spectral resolution. 



7. Results and Conclusions 

We converge atmospheres with effective temperatures T e g = 10 5 ' 5 ,10 6 and 10 6 ' 5 K, with 
vertical (0# = 0) mag netic fields of 10 12 - 5 and 10 13 Gauss. To examine the dependences on 
magnetic geometry, we also converge atmospheres with @b = 45° and 90° at 10 6 K for both field 
strengths. 

Figure 2 gives the temperature structure computed for the Qb = atmospheres. Open 
circles show the locations of the mean radiative optical depths 1 and 10, labeled with the 
average ionization level {Z c r ) at these points. The temperature structure may be compared with 




Fig. 2.— Left: B = 10 12 5 G, Q B = atmosphere structures (p-T plane) for log(T eff ) = 5.5, 6.0, 6.5 
superimposed on curves for Ep/kT = 10~ 4 , 10~ 3 and T = 1, 3 and 10 curves. Positions of r rac j = 1 
and 10 (see §^) are marked with open circles. Labels by the points give the mean ionization level 
at these depths. Right: Same for B = 10 13 G. 



background curves showing V and Ep/kT; the depth zones dominating the spectrum lie in the 
gaseous non-degenerate regime. 

Figure 3 shows the emergent spectra for the normal field atmospheres — these can be 
compared with the corresponding blackbodies. As in Romani (1987), RR96 and ZPS96 we see 
that the rich absorption spectrum of the heavy element gives an overall opacity trend much flatter 
than the steep Kramer's law fall-off of H and He. This assures that the spectra are globally 
closer to a blackbody shape than the light element atmospheres. In fact, the bound-free edges 
at ~ 0.3 — l.OkeV provide a large increase in opacity that gives Fe spectral deficits above the 
Wien peak. Qualitatively, this is similar to the effect of L-edge absorptions in non-magnetic iron 
(RR96). Looking at the opacities and at the absorption features in the spectra, one can see that 
the magnetic Fe transitions are divided, roughly, into two energy ranges: Az^ = transitions 
below ~ lOOeV (for our lower field) and all other transitions, mostly at energies above 300eV. 
Interestingly, for typical neutron star magnetic fields these transitions are at energies roughly 
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comparable to the M- and L-shell transitions in non-magnetic Fe. The resulting spectral shapes 
therefore have some similarities to non-magnetic iron. Note that we do not expect the detailed 
line positions in our atmospheres to be accurate; the general appearance of these emergent spectra 
should however be similar to more accurate results. Although more complete atomic data should 
give an even richer line spectrum, we do find that certain species dominate at depths where the 
lines are forming. Thus high resolution X-ray spectra can in principal extract useful information 
about the Fe atmosphere structure. In particular, the rather strong line broadening, most visible 
at low E, does not prevent discrete lines from being visible at higher energies. At low energies 
lines are also visible, but equivalent widths are questionable because of uncertainty in the detailed 
temperature structure in the very low r surface layers. 




-2 -1 0-2-1 
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Fig. 3.— Left: Emergent Fe spectra for B = 10 12 5 G and Q B = at log(T eff ) = 5.5, 6.0 and 
6.5, showing the two modes (thin solid lines), the total flux (heavy solid line) and corresponding 
black-body spectrum (dashed). Right: Emergent spectra for B = 10 13 G. 

In Figure 3 we also show the contributions of the two normal modes to the emergent flux at 
each temperature. While we do not discuss it in detail here, it is clear that the thermal emission 
will show a strong energy dependent polarization. This polarization also depends on the magnetic 
field geometry; thermal polarization from neutron star surfaces should be pulsed at the star's spin 
period. 
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In practice, we expect quite substantial variation of the emergent spectrum over the surface. 
If due to a coherent dipole, the variation of magnetic field strength and orientation will affect the 
bulk flow of thermal flux from the interior to the star surface. For example Shibanov and Yakovlev 
(1996) indicate an effective temperature decrease as large as ~ 7 from the magnetic pole to the 
magnetic equator for a 10 13 G, (T c g) = 10 6 K neutron star. Such temperature dependence implies 
large thermal emissivity variation with viewing angle. For a given T e g we find that the broadband 
spectral dependence on 6s is quite modest (see also Pavlov, et al. 1995). However, for a surface 
with varying T c g, including the spectral dependence imposed by an atmosphere (e.g. as modeled 
for Geminga by Meyer, Pavlov and Meszaros 1994) can give a complex, energy dependent pulse 
behavior. In fact, spectral features are quite dependent on the magnetic field strength (Figure 3). 
At a minimum, even if the large scale dipole dominates the field at the surface, there will be a 
factor of two difference in B between the magnetic equator and the magnetic poles. Higher order 
multipole contributions are also possible; the large field variations will blur any spectral features 
significantly. 

Finally in addition to the magnetic field variation and effective temperature variation, 
there will be a significant modulation of the neutron star thermal spectrum with changes of the 
observer's viewing angle to an emitting region on the stellar surface. ZPS96 have emphasized the 
'limb darkening' decrease in the emitted flux with increasing [/, for non-magnetic atmospheres. 
In magnetized atmospheres the \i dependence can be extremely complicated in general. One 
basic effect is visible below E 1 = O.lkeV in the normal field case of Figure 4; here one of the 
two normal modes has a significant z polarization component, which experiences a much higher 
continuum opacity unsuppressed by the field. The opacity increase acts in concert with the usual 
limb darkening to diminish the flux at small [i. For tangential Qb = 90°, on the other hand, when 
H ~ in the z-n plane both normal modes are primarily (+, — ), resulting in strongly suppressed 
continuum opacity. The extra transparency counteracts the limb darkening to some extent; the 
emissivity varies little below O.lkeV in the right panel of Figure 4. Thus a simple dipole field 
neutron star would have a relatively bright limb at low E 1 when viewed along the magnetic axis. 

However, other effects can complicate this simplest picture. First, the lines and line wings 
can couple selectively to a single polarization component, and block out the modes transparent to 
the continuum. Secondly, when fi ~ 1, the normal modes are circularly polarized and similar to 
one another, both containing small amounts of the z component; while for [i <^ 0.7 the modes are 
linearly polarized, and the mode perpendicular to the B-k plane is almost completely z-free. Thus 
when the z opacity excess is too extreme, the z-free mode can have larger intensity at small \x 
than both modes together at \i ~ 1. This effect is visible above 0.2 keV in Fig. 4, where strong 
line wings are coupled to the z component. In some cases even the continuum-controlled thermal 
peaks are dominated by this effect; however, any small change in normal mode composition 
can modify this behavior. Clearly models of purely thermal pulsed emission from neutron stars 
can be complex, with effective temperature, line feature position and limb darkening behavior 
controlled by the local magnetic geometry. These effects may be somewhat blurred by gravitational 
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Fig. 4. — Left: Emergent Fe atmosphere emissivity as a function of viewing angle for B = 10 125 G, 
T = 10 5 - 75 K, Q B = 0. Total fluxes are shown at /x = 0.1,0.6 (lighter curve) and 1.0. Note that 
the limb darkening is a strong function of energy. Right: Viewing angle dependence for 6# = tt/2. 
The limb darkening, spectral features and polarization properties are a strong function of magnetic 
geometry. 

defocusing of the emergent radiation (Page 1995), but substantial variation in the atmosphere 
features and overall emissivity are expected as the star rotates. 

Our models suggest that if pulsars possess iron atmospheres, upcoming high throughput, 
moderate resolution spectroscopy missions such as XMM and AXAF will have a wealth of 
spectral features for diagnosis of atmosphere temperature structure and magnetic geometry. 
Data providing such rich spectra will certainly motivate more complete magnetic Fe models. 
Improvement of the atomic data sets used to generate the opacities will be essential for such 
models. Even exhaustive lists of configurations and lines, though, will have to be coupled with 
a more sophisticated treatment of the non-ideal equation of state, along the lines of an activity 
expansion of the grand canonical ensemble (Rogers 1986), or an occupation probability formalism 
(Hummer & Mihalas 1988). Both these approaches, which have been used for non-magnetic 
plasma opacity calculations by the OPAL and OP collaborations respectively, avoid the somewhat 
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ad hoc cutoff procedures necessary in using equation (||). A fully consistent opacity/EOS system 
will also need to include the effects of plasma non-ideality on atomic (b-f and b-b) opacity. This 
will require calculation of states and energies of non-isolated atoms, desirable but prohibitively 
difficult at present. For the line opacities, the dipole selection rules and simple coupling of each 
line to a single polarization component are known to be a first approximation only, subject to 
complication by all forms of broadening. Presumably more detailed atomic modeling will address 
these topics, and produce more precise oscillator strengths as well. Finally, more realistic magnetic 
geometries and inclusion of (weak) mode coupling through scattering will have some effect on the 
emergent flux. Nonetheless, we expect our magnetic models have captured the qualitative behavior 
of heavy element atmospheres. In simplest terms this is that the wealth of spectral features bring 
the overall spectrum shape closer to a blackbody than for light element surfaces. 



Blackbody Fits to Model Data 




Log(T eH ) 

Fig. 5. — Inferred T e fj when fitting a blackbody to a simulated 3000 count ROSAT PSPC spectrum 
produced by a redshifted (z = 0.306), absorbed (Njj = 10 20 cm~ 2 ) model atmosphere. Curves show 
the inferred T for magnetic and non- magnetic H and Fe as a function of the true T e g . 

As in RR96, we can illustrate these broad band spectral effects by simulating a modest 
statistic (3000 count) PSPC exposure of a neutron star, whose absorbed spectrum is fitted by 
a blackbody. Comparing the inferred Tbb with the true T e g for various atmosphere models 
illustrates the low resolution spectral effects (Figure 5). As for non-magnetic atmospheres, the 
Fe results are much softer for a given T e g- than H atmospheres. In contrast, existing data suggest 
phase average spectra appreciably harder than the corresponding blackbody. While the clear 
existence of magnetospheric flux contributions and the likely T c g variations over the surface 
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complicate the interpretation, it would seem that these data support the existence of light element 
atmospheres on neutron star surfaces. If substantiated, this conclusion would give important 
clues to the history of the neutron star surface and serious implications for cooling curves and 
the equation of state of matter at high density. The magnetic model Fe atmospheres produced in 
this paper are of sufficient quality to allow selection between H models and Fe models for pulsar 
X-ray surface emission. If H atmospheres are preferred, the sophisticated modeling possible for 
this simple species will allow detailed interpretation of soft X-ray spectra. If however, Fe models 
remain acceptable, further careful study of the behavior of highly magnetized iron under neutron 
star surface conditions will be strongly motivated. 

We thank G. G. Pavlov for supplying the sample magnetic H model spectra and for giving 
a careful and detailed critique of this paper. This work was supported in part by the National 
Sciences and Engineering Research Council of Canada (MR); by NASA grants NAG 5-3101, and 
NAGW-4526 (RWR); and by NAG 5-2868 and, through the Compton fellowship program, by 
NASA grant NAG 5-2687 (MCM). 
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